R/funcion polinomios 150419.R

In nxarteagag/rddsw: Regression discontinuity with and without sample weights

#' @title RDD poliynomial local regression
#' @description no parametric regression discontinuity using polynomial local regression estimator with and without sample weights
#' @param D variable with: 1 for treatment group and 0 for control group
#' @param Z continuous variable
#' @param Y outcome variable
#' @param c cutpoint
#' @param type sharp or fuzzy
#' @param TreatLeft T if the treatment is applied for below of the cutpoint
#' @param W sample weights 1/pik, NULL if there aren't sample weights
#' @param svydesign_ design sample using survey library
#' @param vcov_ variances and covariances type: HCO, HC1, HC2 y HC3
#' @param kernel_ kernel type: triangular, epanechnikov, tricube, uniform, gaussian y cosine
#' @param h bandwidth
#' @param q order polynomial
#' @return Impact estimation, e.e. and t test
#' @examples
#'
RDpol <- function(D = muestra1$D, Z = muestra1$x, Y = muestra1$y,c =0,type="sharp",TreatLeft=T,weights=NULL,svydesign_=muestra,vcov_,kernel_="uniforme",h=2,q=2){
  library(survey)
  if(!is.null(weights)){
    muestra        <- data.frame(D,Z,Y,weights)
  }else{
    muestra        <- data.frame(D,Z,Y)
  }
  muestra$intercepto  <- 1
  muestra$Z_C         <- muestra$Z-c
  muestra$ZD          <- muestra$Z*muestra$D
  if(kernel_!="gaussian"){
    muestra <- muestra[which(muestra$Z>=c-h & muestra$Z<=c+h),]
    if(!is.null(weights)){
      svydesign_ <- subset(svydesign_,Z>=c-h & Z<=c+h)
    }
  }
  Y        <- muestra$Y

  if(!is.null(weights)){
    W <- kernel_sRD(kernel_=kernel_,weights=muestra$weights,h=h,X=muestra$Z,c=c)
  }else{
    W <- kernel_sRD(kernel_=kernel_,weights=NULL,h=h,X=muestra$Z,c=c)
  }
  if(type!="sharp"){
    muestra$T_          <- ifelse(muestra$Z<c,1,0)
    if(TreatLeft==F){
      muestra$T_          <- 1-muestra$T_
    }
    muestra$ZT_          <- muestra$Z*muestra$T_
  }
  if(type=="sharp"){
    Xm1                 <- as.matrix(muestra[,c("intercepto","Z_C","D","ZD")])
    for(i in 2:q){
      Xm1 <- cbind(Xm1,Xm1[,2]^i,Xm1[,4]^i)
      colnames(Xm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
    }
    impacto.1_2         <- solve(t(Xm1)%*%W%*%Xm1)%*%(t(Xm1)%*%W%*%Y)
    beta                <- impacto.1_2[3,]
  }else{if(type=="fuzzy"){
    Zm1            <- as.matrix(muestra[,c("intercepto","Z_C","T_","ZT_")])
    for(i in 2:q){
      Zm1 <- cbind(Zm1,Zm1[,2]^i,Zm1[,4]^i)
      colnames(Zm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
    }
    Xm1            <- as.matrix(muestra[,c("intercepto","Z_C","D","ZT_")])
    for(i in 2:q){
      Xm1 <- cbind(Xm1,Xm1[,2]^i,Xm1[,4]^i)
      colnames(Xm1)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
    }

    beta1stg       <- solve(t(Zm1)%*%W%*%Zm1)%*%(t(Zm1)%*%W%*%muestra$D)
    muestra$Xest   <- Zm1%*%beta1stg
    Xest           <- as.matrix(muestra[,c("intercepto","Z_C","Xest","ZT_")])
    for(i in 2:q){
      Xest <- cbind(Xest,Xest[,2]^i,Xest[,4]^i)
      colnames(Xest)[(3+2*(i-1)):(4+(2*(i-1)))] <- c(paste0("Z_C_",q),paste0("ZD_",q))
    }
    impacto.1_2    <- solve(t(Xest)%*%W%*%Xest)%*%(t(Xest)%*%W%*%Y)
    beta           <- impacto.1_2[3,]
  }
  }

  # Varianza sin pesos
  Yest                <- Xm1%*%impacto.1_2
  sigma2              <- sum((Y - Yest)^2)/(dim(muestra)[1]-length(impacto.1_2))
  M_                  <- (Y - Yest)^2
  M                   <- matrix(0,length(M_),length(M_))
  diag(M)             <- M_
  if(type=="sharp"){
    if(is.null(weights)){
      if(vcov_=="HC0"){
        var              <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
      }else{if(vcov_=="HC1"){
        mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
        diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
        var       <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
      }else{if(vcov_=="HC2"){
        mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
        diag(mat) <- 1/((1-diag(Xm1%*%solve(t(Xm1)%*%Xm1)%*%t(Xm1))))
        var       <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
      }else{if(vcov_=="HC3"){
        mat <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])
        diag(mat) <- 1/((1-diag(Xm1%*%solve(t(Xm1)%*%Xm1)%*%t(Xm1)))^2)
        var       <- solve(t(Xm1)%*%W%*%Xm1)%*%t(Xm1)%*%W%*%mat%*%M%*%t(W)%*%Xm1%*%solve(t(Xm1)%*%W%*%Xm1)
      }}}}
    }else{
      var1            <- solve(t(Xm1)%*%W%*%Xm1) # revisar si dejo o pongo sigma2
      muestra$ek2     <- muestra$Y-Xm1%*%impacto.1_2
      diag(M)         <- muestra$ek2
      mat             <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])

      if(vcov_=="HC0"){
        diag(mat) <- 1
      }else{if(vcov_=="HC1"){
        diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
      }}

      svydesign_$variables$uk <- mat%*%M%*%t(W)%*%Xm1
      varmuestra           <- vcov(svytotal(~uk,svydesign_))
      var                  <- (var1)%*%varmuestra%*%(var1)
    }
  }else{if(type=="fuzzy"){
    if(is.null(weights)){
      if(vcov_=="HC0"){
        var       <- solve(t(Xest)%*%W%*%Xest)%*%t(Xest)%*%W%*%M%*%t(W)%*%Xest%*%solve(t(Xest)%*%W%*%Xest)
      }else{if(vcov_=="HC1"){
        mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
        diag(mat) <- (dim(Xest)[1]/(dim(Xest)[1]-length(impacto.1_2)))
        var       <- solve(t(Xest)%*%W%*%Xest)%*%t(Xest)%*%W%*%mat%*%M%*%t(W)%*%Xest%*%solve(t(Xest)%*%W%*%Xest)
        #var       <- (dim(Xest)[1]/(dim(Xest)[1]-2))*var
      }else{if(vcov_=="HC2"){
        mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
        diag(mat) <- 1/((1-diag(Xest%*%solve(t(Xest)%*%Xest)%*%t(Xest))))
        var       <- solve(t(Xest)%*%Xest)%*%t(Xest)%*%mat%*%M%*%Xest%*%solve(t(Xest)%*%Xest)
      }else{if(vcov_=="HC3"){
        mat <- matrix(0,dim(Xest)[1],dim(Xest)[1])
        diag(mat) <- 1/((1-diag(Xest%*%solve(t(Xest)%*%Xm1)%*%t(Xest)))^2)
        var       <- solve(t(Xest)%*%Xest)%*%t(Xest)%*%mat%*%M%*%Xest%*%solve(t(Xest)%*%Xest)
      }}}}
    }else{
      var1            <- solve(t(Xest)%*%W%*%Xest) # revisar si dejo o pongo sigma2
      muestra$ek2     <- muestra$Y-Xm1%*%impacto.1_2
      diag(M)         <- muestra$ek2
      mat             <- matrix(0,dim(Xm1)[1],dim(Xm1)[1])

      if(vcov_=="HC0"){
        diag(mat) <- 1
      }else{if(vcov_=="HC1"){
        diag(mat) <- (dim(Xm1)[1]/(dim(Xm1)[1]-length(impacto.1_2)))
      }}
      svydesign_$variables$uk <- mat%*%M%*%t(W)%*%Xm1
      varmuestra              <- vcov(svytotal(~uk,svydesign_))
      var                     <- (var1)%*%varmuestra%*%(var1)
    }
  }}

  seBeta               <- sqrt(var[3,3])
  significativo        <- ifelse(abs(beta[i]/(seBeta[i]))>=qt(0.975,dim(muestra)[1]-dim(Xm1)[2]),"Significativo","No Significativo")
  impact               <- data.frame(round(beta,3),round(seBeta,5),pvalue=(1-pt(abs(beta/(seBeta)),df = dim(muestra)[1]-length(impacto.1_2)))*2)
  colnames(impact)     <- c("Impacto","S.E.","P-Value")
  impact$significativo <- ifelse(impact$`P-Value`<0.01,"***",ifelse(impact$`P-Value`<0.05,"**",ifelse(impact$`P-Value`<0.1,"*","")))
  colnames(impact)[4]  <- ""
  print(impact)
}